library(readxl)
library(ggplot2)
library(grid)
library(gridExtra)
library(png)
library(sf)
library(spatialreg)
library(spdep)
library(ncdf4)
library(RColorBrewer)
library(cowplot)
library(colorspace)
library(captioner)
fig_nums <- captioner()
table_nums <- captioner(prefix = "Table")
memory.limit(size = 160000)
## [1] Inf
The TOC concentration from the 1000-lakes survey in 2019 was provided by NIVA (NIVA2020?) and is available on GitHub (Crapart n.d.). The lakes samples in 2019 were matched with the lakes sampled in 1995 based on their sampling coordinates, then the data was extracted in the same way as for the model with 1995-data.
The TOC data and sampling points were provided by Heleen de Wit.
niva <- read_xlsx("NIVA_test/niva_selection.xlsx")
The lakes from the 1000-lakes survey were matched with the lakes from the Northern European Lake Survey based on the sampling point. The correspondance between 1995 and 2019 TOC-values was checked in a second step.
fennoscandia.33 <- readRDS("fennoscandia.33.rds") %>% dplyr::select(c("ebint","nation","TOC","longitude","latitude","geom"))
niva <- read_xlsx("NIVA_test/niva_selection.xlsx") %>% dplyr::select(c("station_id","longitude","latitude","toc_2019","toc_1995"))
niva.sf <- st_as_sf(niva,coords=c("longitude","latitude"))
st_crs(niva.sf) <- "+proj=longlat +datum=WGS84 +no_defs"
niva.sf <- niva.sf %>% st_transform(crs(fennoscandia.33))
fennoscandia.point <- fennoscandia.33 %>% st_drop_geometry() %>% st_as_sf(coords = c("longitude","latitude"), remove = F)
st_crs(fennoscandia.point) <- "+proj=longlat +datum=WGS84 +no_defs"
fennoscandia.sf <- st_transform(fennoscandia.point,crs(fennoscandia.33))
niva.near <- st_join(st_buffer(niva.sf, dist = 1000), st_buffer(fennoscandia.sf, 1000), join = st_intersects)
niva.unique <- niva.near %>% dplyr::filter(toc_1995 == TOC)
niva.poly <- merge(dplyr::select(fennoscandia.33,c("ebint","geom")), st_drop_geometry(niva.unique), by = "ebint")
saveRDS(niva.poly,"NIVA_test/niva.poly.rds")
niva.poly <- readRDS("NIVA_test/niva.poly.rds")
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(st_crs(niva.poly))
map <- ggplot()+
geom_sf(data = nor, fill = "white", size = 0.2)+geom_sf(data = niva.poly, aes(fill = toc_2019, col = toc_2019))+
scale_color_distiller(type = "div", palette = "RdYlBu",name = "TOC sampled 2019 (mg/L)", aesthetics = c("colour","fill"))+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = map, filename = "NIVA_test/map.TOC.2019.png", dpi = "retina", width = 10, height = 12, units = "cm")
compare <- ggplot(niva.poly)+geom_point(aes(x=TOC,y=toc_2019, col = latitude), size = 1)+
scale_color_continuous_divergingx(palette = "RdYlBu", name = "Latitude", mid = 65)+
labs(x = "TOC concentration in lakes,\n 1995 (mg/L)", y = "TOC concentration in lakes, 2019 (mg/L)")+
geom_abline(slope = 1, intercept = 0, col = "black")+
annotate(geom = "text", label = "1:1", x = 12, y = 11, size = 3)+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = compare, filename = "NIVA_test/compare.toc.png", dpi = "retina", width = 10, height = 12, units = "cm")
listgrob <- list("NIVA_test/map.TOC.2019.png","NIVA_test/compare.toc.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob,ncol = 2, nrow = 1, labels = c("a)","b)"), label_size = 20)
The NDVI was extracted for the year 2015 from the GIMMS dataset (The National Center for Atmospheric Research 2018), as data from this database was not available for later dates. Other databases like Copernicus could have provided this data, but GIMMS was used for the original model and we assumed that the changes in mean NDVI are minimal in 3 years (as NDVI is taken for y-1). For details about the extraction of NDVI in each catchment, see Supplementary 1.
# http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/
library(gimms)
ndvi.2015 <- downloadGimms(x= 2015,y=2015,dsn = "NDVI")
ndvi.max.2015 <- monthlyComposite(ndvi.2015,monthlyIndices(ndvi.2015))
saveRDS(ndvi.max.2015,"NIVA_test/ndvi.max.2015.rds")
ndvi.max.2015 <- readRDS("NIVA_test/ndvi.max.2015.rds")
niva.poly <- readRDS("NIVA_test/niva.poly.rds") %>% dplyr::select(c("ebint","geometry"))
# Makes the raster smaller, ensuring faster processing of the data
summer.scandinavia.2015 <- raster::crop(ndvi.max.2015,c(0,35,55,73))
# Re-introduces NA
summer.scan <- reclassify(summer.scandinavia.2015, cbind(-Inf, 0, NA), right=FALSE)
# Stack bands for summer
summer.mean.2015 <- raster::stack(summer.scan[[6]],summer.scan[[7]],summer.scan[[8]]) %>% mean()
# extract data
summer.ndvi.2015 <- raster::extract(summer.mean.2015,niva.poly, fun = mean, df = T, sp = T, na.rm = T)
names(summer.ndvi.2015) <- c("ebint","ndvi")
summer.ndvi.2015$NDVI <- floor(summer.ndvi.2015$ndvi/10)/1000
# save df
saveRDS(summer.ndvi.2015, "NIVA_test/100lakes.summer.ndvi.2015.rds")
write.csv(summer.ndvi.2015,"NIVA_test/100lakes.summer.ndvi.2015.csv")
knitr::include_graphics("NIVA_test/map.ndvi.NIVA.png")
The runoff data was downloaded from the CORDEX database (CORDEX 2021). We used the average over the interval 1985-2015 (last available year in historical models).
niva.poly <- readRDS("NIVA_test/niva.poly.rds") %>% dplyr::select(c("ebint","geometry"))
runoff.raster <- "CORDEX/mrros_Lmon_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_185001-201412.nc"
runoff.nc <- nc_open(runoff.raster)
runoff.stack <- raster::stack(runoff.raster, bands = c(1621:1980)) # runoff from 1985 to 2015
runoff.mean <- runoff.stack %>% mean()
runoff.niva <- raster::extract(runoff.mean, niva.poly, fun = mean, df = T, sp = T, na.rm = T)
names(runoff.niva) <- c("ebint","Runoff")
saveRDS(runoff.niva,"NIVA_test/runoff.niva.rds")
ggplot(st_as_sf(runoff.niva))+geom_sf(aes(fill=Runoff,col=Runoff))
knitr::include_graphics("NIVA_test/map.runoff.NIVA.png")
The Land Cover data was downloaded from https://land.copernicus.eu/pan-european/corine-land-cover/clc2018. We used the 2018 version (last one available).
The land cover data was downloaded as a raster (.tiff file), which included the legend recording the code for each land use. The raster and legend were assembled in QGIS before being imported as a shapefile in R. The codes (corresponding to the line in the extracted dataframe) of the categories of interest were:
Arable land:
Bogs:
Forest:
Bare:
31: Bare rocks 32: Sparsely vegetated areas 33: Burnt areas
niva.poly <- readRDS("NIVA_test/niva.poly.rds") %>% dplyr::select(c("ebint","geometry"))
corine.2018 <- raster::raster("NIVA_test/CLC/U2018_CLC2018_V2020_20u1.tif")
niva.corine <- niva.poly %>% st_transform(crs(corine.2018))
niva.list <- split(niva.corine,seq(nrow(niva.corine)))
extract_corine <- function(x){
clc <- extract(corine.2018,x)
saveRDS(clc,file = paste("NIVA_test/CLC/clc",x$ebint,"rds",sep="."))
}
lapply(niva.list,extract_corine)
clc.test <- readRDS("NIVA_test/CLC/clc.16711649.rds")
clc.files <- list.files(path = "NIVA_test/CLC/", pattern = "clc.*.rds", full.names = T)
list.clc <- sapply(clc.files,readRDS)
list.clc.niva <- sapply(list.clc, as.integer)
names.list <- str_extract(clc.files,"\\d+")
names(list.clc.niva) <- names.list
saveRDS(list.clc.niva,"NIVA_test/CLC/list.clc.niva.rds")
list.clc.niva <- readRDS("NIVA_test/CLC/list.clc.niva.rds")
clc.tab.niva <- sapply(list.clc.niva, function(x) tabulate(x,45)) %>% t()
clc.tab.area.niva <- clc.tab.niva*prod(res(corine.2018))
catchment.area <- rowSums(clc.tab.area.niva)
clc.tab.prop.niva <- sweep(clc.tab.area.niva,1,catchment.area, FUN = "/")
saveRDS(clc.tab.prop.niva,"NIVA_test/CLC/clc.tab.prop.niva.rds")
clc.tab.prop.niva <- readRDS("NIVA_test/CLC/clc.tab.prop.niva.rds") %>% as.data.frame()
bogs <- clc.tab.prop.niva[,36] %>% as.data.frame() %>% setNames("Bog") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(bogs,"NIVA_test/CLC/bogs.niva.rds")
arable <- rowSums(clc.tab.prop.niva[,c(12,16,18,19,20)]) %>% as.data.frame() %>% setNames("Arable") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(arable,"NIVA_test/CLC/arable.niva.rds")
forest <- rowSums(clc.tab.prop.niva[,c(23,24,25)]) %>% as.data.frame() %>% setNames("Forest") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(forest,"NIVA_test/CLC/forest.niva.rds")
bare <- rowSums(clc.tab.prop.niva[,c(31,32,33)]) %>% as.data.frame() %>% setNames("Bare") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(bare,"NIVA_test/CLC/bare.niva.rds")
listgrob <- list("NIVA_test/map.Bog.NIVA.png","NIVA_test/map.arable.NIVA.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 1, align = "v", greedy = T, labels = c("a)","b)"), label_size = 20)
The nitrogen and sulfur deposition data were extracted from the EMEP database (https://emep.int/mscw/mscw_moddata.html#Comp). The data from 2019 was used.
N-deposition (NOx and NH3) is the sum of: * Dry deposition of oxidized nitrogen per m2 grid DDEP_OXN_m2Grid, in mg/m2 * Wet deposition of oxidized nitrogen WDEP_OXN, in mg/m2 * Dry deposition of oxidized nitrogen per m2 grid DDEP_RDN_m2Grid, in mg/m2 * Wet deposition of reduced nitrogen WDEP_RDN, in mg/m2
EMEP_file <- "NIVA_test/Ndep/EMEP01_rv4.42_year.2019met_2019emis.nc"
niva.poly <- readRDS("NIVA_test/niva.poly.rds") %>% dplyr::select(c("ebint","geometry"))
woxn <- raster(EMEP_file, varname = "WDEP_OXN")
doxn <- raster(EMEP_file, varname = "DDEP_OXN_m2Grid")
wrdn <- raster(EMEP_file, varname = "WDEP_RDN")
drdn <- raster(EMEP_file, varname = "DDEP_RDN_m2Grid")
woxn_df <- extract(woxn,niva.poly, sp = T, df = T, na.rm = T, fun = mean)
names(woxn_df) <- c("ebint","woxn")
woxn_u <- woxn_df %>% st_as_sf() %>% st_drop_geometry() %>% dplyr::distinct()
doxn_df <- extract(doxn,niva.poly, sp = T, df = T, na.rm = T, fun = mean)
names(doxn_df) <- c("ebint","doxn")
doxn_u <- doxn_df %>% st_as_sf() %>% st_drop_geometry() %>% dplyr::distinct()
wrdn_df <- extract(wrdn,niva.poly, sp = T, df = T, na.rm = T, fun = mean)
names(wrdn_df) <- c("ebint","wrdn")
wrdn_u <- wrdn_df %>% st_as_sf() %>% st_drop_geometry() %>% dplyr::distinct()
drdn_df <- extract(drdn,niva.poly, sp = T, df = T, na.rm = T, fun = mean)
names(drdn_df) <- c("ebint","drdn")
drdn_u <- drdn_df %>% st_as_sf() %>% st_drop_geometry() %>% dplyr::distinct()
ndep.df <- merge(woxn_u,doxn_u, by = "ebint") %>% merge(wrdn_u, by = "ebint") %>% merge(drdn_u, by = "ebint")
ndep.df$TNdep <- rowSums(ndep.df[,c(2:5)])
saveRDS(ndep.df,"NIVA_test/Ndep/ndep.df.niva.rds")
knitr::include_graphics("NIVA_test/map.TNdep.NIVA.png")
The dataset from NIVA and the NDVI values were merged and the variables were given names matching the Fennoscandian dataset. The projection of the polygons was converted to the projection used in the present project (EU33).
ndvi.niva <- readRDS("NIVA_test/100lakes.summer.ndvi.2015.rds") %>% st_as_sf() %>% st_drop_geometry() %>% dplyr::distinct()
runoff.niva <- readRDS("NIVA_test/runoff.niva.rds") %>% st_as_sf() %>% st_drop_geometry() %>% dplyr::distinct()
bogs.niva <- readRDS("NIVA_test/CLC/bogs.niva.rds")
arable.niva <- readRDS("NIVA_test/CLC/arable.niva.rds")
ndep.df.niva <- readRDS("NIVA_test/Ndep/ndep.df.niva.rds")
niva.poly <- readRDS("NIVA_test/niva.poly.rds")
niva.data <- niva.poly %>% merge(ndvi.niva, by = "ebint") %>% merge(runoff.niva, by = "ebint") %>% merge(bogs.niva, by = "ebint") %>% merge(arable.niva, by = "ebint") %>% merge(ndep.df.niva, by = "ebint") %>% dplyr::distinct()
niva.data <- niva.data[-which(is.na(niva.data$toc_2019) == T),]
niva.data <- niva.data[-which(is.na(niva.data$NDVI) == T),]
niva.data$Runoff <- niva.data$Runoff *365*24*60*60 # from kg/m2/s to kg/m2/year or mm/y
niva.data$logRunoff <- log(niva.data$Runoff)
saveRDS(niva.data,"NIVA_test/niva.data.rds")
niva.data <- readRDS("NIVA_test/niva.data.rds")
m <- c("TOC","NDVI","Runoff","Bog","Arable","TNdep")
u <- c("TOC concentration,\nmgC/L", "NDVI on scale 0 to 1", "Mean Runoff, mm/y","Proportion of bogs","Proportion of\narable land","Total nitrogen\ndeposition, mg/m2")
dir <- c(T,F,F,T,T,F)
midpoint <- (c(median(niva.data$TOC),median(niva.data$NDVI),median(niva.data$Runoff),0.1,0.1,median(niva.data$TNdep)))
marg <- 0
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(crs(niva.data))
swe <- st_read("Country_shapefile/sweden.shp") %>% st_transform(crs(niva.data))
fin <- st_read("Country_shapefile/finland.shp") %>% st_transform(crs(niva.data))
for(i in m){
mysf <- niva.data
midcol <- median(mysf[[i]])
map <- ggplot()+geom_sf(data = nor, fill = "white", size = 0.2)+
geom_sf(data = mysf, aes(fill = .data[[i]], col = .data[[i]]))+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name= paste(u[grep(i,m)],"\n(",i,")",sep=""), rev = dir[grep(i,m)], mid = midpoint[grep(i,m)]) +
theme_minimal(base_size = 8)+
theme(legend.position = "bottom", plot.margin = unit(c(0,marg,0,marg), units = "cm"))
filename_i <- paste("NIVA_test/map",i,"NIVA","png",sep = ".")
ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
}
In order to predict the TOC in 2019 using the SELM, a neighbor matrix with the NIVA-lakes is first computed.
niva.data <- readRDS("NIVA_test/niva.data.rds")
niva.kmat <- st_centroid(niva.data, of_largest_polygon = T) %>% knearneigh(k = 100) %>% knn2nb() %>% nb2listw()
saveRDS(niva.kmat,"NIVA_test/niva.kmat.rds")
The SELM fitted in Supplementary 1, with NDVI, Bog, Arable, Runoff and TNdep as predictors, is then applied to the 2019-dataset to predict the TOC concentration in lakes. The difference between model results and observations is computed as \(ln([TOC]_{obs})-ln([TOC]_{fit})\).
sem.model <- readRDS("sem.fennoscandia.5.rds")
niva.data <- readRDS("NIVA_test/niva.data.rds")
niva.kmat <- readRDS("NIVA_test/niva.kmat.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
niva_fit <- predict(sem.model, newdata = niva.data, niva.kmat)
niva_predict <- cbind(niva.data, niva_fit)
niva_predict$diffTOC <- niva_predict$toc_2019 - exp(niva_predict$fit)
niva_predict$diffLogTOC <- log(niva_predict$toc_2019) - niva_predict$fit
saveRDS(niva_predict,"NIVA_test/niva_predict.rds")
niva_predict <- readRDS("NIVA_test/niva_predict.rds")
ggplot(niva_predict) + geom_point(aes(y = log(toc_2019), x = fit, col = latitude), size = 1)+
scale_color_continuous_divergingx(palette = "RdYlBu", aesthetics = c("col"), mid = 60, name = "Latitude")+
labs(y = "logTOC 2019", x = "Predicted logTOC")+
annotate(geom = "label", label = paste("r = ", round(cor(log(niva_predict$toc_2019),niva_predict$fit),2)), x = 2, y = -1.5, size = 3)+
geom_abline(slope = 1, intercept = 0, col = "black")+
annotate(geom = "text", label = "1:1", x = 2.5, y = 2.2, size = 3)+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(filename = "NIVA_test/niva-test-result.png", dpi = "retina", width = 10, height = 12, units = "cm")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
niva_predict <- readRDS("NIVA_test/niva_predict.rds")
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(st_crs(wr.sf.95))
mypal <- brewer.pal(9, name = "RdYlBu")
lims.toc <- c(min(c(log(niva_predict$toc_2019), niva_predict$fit)), max((c(log(niva_predict$toc_2019), niva_predict$fit))))
diff.niva <- ggplot() + geom_sf(data = nor, fill = "white", col = "black", size = 0.1)+
geom_sf(data = niva_predict, aes(fill = diffLogTOC, col = diffLogTOC)) +
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="Difference between\nobserved and predicted logTOC\nwith NDVI, Bog, Arable, logRunoff and TNdep", mid = 0, limits = lims.toc, rev = F) +
labs(x = "", y = "") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom")
ggsave(plot = diff.niva, filename = "NIVA_test/diff.niva.png", dpi = "retina", width = 10, height = 12, units = "cm")
diff.niva.sub <- ggplot() + geom_sf(data = nor, fill = "white", col = "black", size = 0.1)+
geom_sf(data = niva_predict, aes(fill = diffLogTOC, col = diffLogTOC)) +
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="Difference between\nobserved and predicted logTOC\nwith NDVI, Bog, Arable, logRunoff and TNdep", mid = 0, limits = lims.toc, rev = F) +
labs(x = "", y = "", subtitle = "c)") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom")
ggsave(plot = diff.niva.sub, filename = "NIVA_test/diff.niva.sub.png", dpi = "retina", width = 10, height = 12, units = "cm")
listgrob <- list("NIVA_test/niva-test-result.png","NIVA_test/diff.niva.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob,ncol = 2, nrow = 1, labels = c("a)","b)"), label_size = 20)
For comparison, we try to model current TOC concentration based on a linear model. The results are shown below. The difference between model results and observations is computed as \(ln([TOC]_{obs})-ln([TOC]_{fit})\).
lm.model <- readRDS("lm.fennoscandia.5.rds")
niva_final <- readRDS("NIVA_test/niva_final.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
niva_fit <- predict(lm.model, newdata = niva_final)
niva_predict <- cbind(niva_final, niva_fit)
niva_predict$diffTOC <- niva_predict$toc_2019 - exp(niva_predict$niva_fit)
niva_predict$diffLogTOC <- log(niva_predict$toc_2019) - niva_predict$niva_fit
saveRDS(niva_predict,"NIVA_test/niva_predict_lm.rds")
niva_predict <- readRDS("NIVA_test/niva_predict_lm.rds")
ggplot(niva_predict) + geom_point(aes(x = log(toc_2019), y = niva_fit, col = latitud))+
labs(x = "log(TOC 2019)", y = "Predicted logTOC")+
scale_color_continuous_divergingx(palette = "RdYlBu", name = "Latitude", mid = 60)+
annotate(geom = "label", label = paste("r = ", round(cor(log(niva_predict$toc_2019),niva_predict$niva_fit),2)), x = 2, y = -0.2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(filename = "NIVA_test/niva-lm-result.png", dpi = "retina", width = 10, height = 12, units = "cm")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
niva_predict <- readRDS("NIVA_test/niva_predict_lm.rds")
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(st_crs(wr.sf.95))
lims.toc <- c(min(c(log(niva_predict$toc_2019), niva_predict$fit)), max((c(log(niva_predict$toc_2019), niva_predict$fit))))
diff.niva <- ggplot() + geom_sf(data = nor, fill = "white", col = "lightgray")+
geom_sf(data = niva_predict, aes(fill = diffLogTOC, col = diffLogTOC)) +
scale_color_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Difference between\nobserved and predicted logTOC\nwith NDVI, Bogs, Arable, logRunoff and TNdep ", aesthetics = c("fill","col"))+
labs(x = "", y = "")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = diff.niva, filename = "NIVA_test/diff.niva.lm.png", dpi = "retina", width = 10, height = 12, units = "cm")
listgrob <- list("NIVA_test/niva-lm-result.png", "NIVA_test/diff.niva.lm.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
grid.arrange(grobs = listgrob, ncol = 2, nrow = 1, labels = c("a)","b)"))
r knit::exit()